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Abstract 

We extend the study of lowest moments, {x) and (x^), of the parton distribution function of the nucleon 
to include those of the sea quarks; this entails a disconnected insertion calculation in lattice QCD. This is 
carried out on a 16^ x 24 quenched lattice with Wilson fcrmion. The quark loops are calculated with Z2 
noise vectors and unbiased subtractions, and multiple nucleon sources are employed to reduce the statistical 
errors. We obtain 5(t signals for {x) for the u,d, and s quarks, but (x^) is consistent with zero within errors. 
We provide results for both the connected and disconnected insertions. The perturbatively renormalized 
{x) for the strange quark at /i = 2 GeV is {x)s+s — 0.027 ± 0.006 which is consistent with the experimental 
result. The ratio of (x) for s vs. u/d in the disconnected insertion with quark loops is calculated to be 
0.88 ± 0.07. This is about twice as large as the phenomenologically fitted -;— ; — ^ ^ from experiments 

{X)u + {X)d 

where u and d include both the connected and disconnected insertion parts. We discuss the source and 
implication of this difference. 

PACS numbers: 11. 15. Ha, 12.28.Gc, ll.30.Rd 



I. INTRODUCTION 



Recently, there has been a good deal of interest in the study of sea quarks, both in theory and 
experiment. In such studies, strange quarks play an important role in observables involving sea 
quarks. For more than a decade, intensive studies have been made in measuring and understanding 
the strangeness contribution to the nucleon spin [1], the electromagnetic form factors [2, 3, 4, 5], the 
strangeness condensate [6], and the parton distribution function in the nucleon [7, 8, 9]. Another 
important aspect of studying strangeness content is to address the issue of the NUTEV anomaly. 
The NUTEV experiment [v + N ^ (i+X) [10], which measures the Weinberg angle or weak mixing 
angle, an important parameter in the Standard Model of particle physics, finds a value which is three 
standard deviations away from the world average value. One suggestion to explain this discrepancy 
is the asymmetry in strange and anti-strange parton distribution [11, 13], a non-perturbative 
effect. Attempts by various theoretical models to calculate the asymmetry give inconsistent results 
(sometimes with completely opposite signs) [13, 14, 15, 16, 17]. Also, phenomenological extractions 
by CTEQ and NUTEV [19, 20] for leading order and next-to-leading order give different results. 
Lattice QCD can assess this asymmetry from first principles in terms of {x^)s-s to address whether 
it is strange or anti-strange which is leading in large x. This information will be helpful to constrain 
and analyze the experimental data. Similarly, the first moment of the strange parton distribution, 
{x)s+s is not well known. It ranges between 0.018 and 0.04 from the fitting of parton distribution 
functions to experiments [8]. 

In the present work, we will study the first and second moments of quark distribution for up, 
down and strange quarks. The first moment provides the measure of the symmetric contribution 
from parton and anti-parton distributions {q + q) and the second moment provides the measure of 
the asymmetry in parton and anti-parton distributions {q — q)- These moments have contributions 
both from connected and disconnected insertions for up and down quarks and only disconnected 
insertion for strange quarks. Since lattice calculations on connected insertions, for the first and 
second moments, have been done before, this work is going to focus mainly on disconnected insertion 
contributions (particularly for strange quarks) to the first and second moments, which has not been 
attempted in lattice QCD. 

This paper is organized as follows. We give the formalism and lattice operators in Sec. II. The 
disconnected insertion calculation is presented in Sec. III. The perturbative renormalization is 
given in Sec. IV. Sec. V presents numerical parameters and error studies of the noise estimate. The 
results for both the disconnected insertions and connected insertions are given in Sec. VI. Finally, 
we offer a conclusion and some discussion in Sec. VII. Some details of the three-point correlation 
functions are given in the Appendices. 

II. FORMALISM 

In deep inelastic scattering [22] , it is useful to consider and analyze the moments of the structure 
function via the operator product expansion (OPE) where, in the limit of distance a; — > or 
equivalently q — > oo, the product of two operators can be expanded in terms of local operators 



where c's are the Wilson coefficients. 

The leading term for such an expansion has the lowest twist, t = 2. For unpolarized structure 
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functions with vector currents J^, the twist-two operators for quarks have the bihnear form 

where ipf denotes the quark field operator for the flavor f, V = \{V — V), and {• • • } stands for 
symmetrization of the indices, /i's. The subtracted trace terms are proportional to gf^^f^j ^ so that 
the operator is traceless on all pairs of indices. 

In the leading twist, the moments of structure functions Fi and F2 can be written as 

/ 

(for n > 2) (3) 

/ 

where is defined through the forward matrix elements 

( P I oWa*! -/^- I p ^ ^ 2 P^i • • • P^" - traces. (4) 

In the parton model, has the interpretation as the (n — l)th moment of the momentum fraction 
carried by the quarks with flavor / at some scale i.e. 

A^^iy.) = [\xx^-' [fix) + (-!)"/»] , (5) 
Jo 

where f{x) is the quark distribution function and f{x) is the anti-quark distribution function for 
the flavor /. We see that the first moment (n = 2) has the symmetric combination of the quark 
and anti-quark distribution and the second moment {n = 3), due to the interference between the 
vector and axial- vector part of the weak interaction current, has the asymmetric combination of 
the quark and anti-quark distribution. Our goal is to compute the first and second moments for 
up, down, and strange quarks. 



2 [ clxx""-^ Fi{x,Q^) 
Jo 

!\xx^-^F2{x,Q^) 
Jo 



A. Lattice Operators 



Lattice calculations are carried out in Euclidean path-integral. Thus, we need to transform the 
twist-two operators from Minkowski space to Euclidean space. Following the convention [23, 24] 

■j^mo _^ _d(e) .^(m),- _^ 

with the 7 matrices defined as 

{l^l.lv} = 2V , 75 = 71727374 , (y^iv = ^ [l^i.lu], (7) 
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we can transform the twist-two operators in Eq. (2) to its Euclidean counterpart by using the 
notation [25] 



{M){n)ui---iyn 



(8) 



where hfj_u = diag (z, i, i, 1) and df^u = diag (— 1, — 1, — 1, i). We have also set a^'^^'^ = a^^^^^ for any 
four vector, a. From now on we will consider Euclidean operators only and drop the superscript, 
E. Since the Euclidean signature is (H — I — I — h), the subtracted trace terms are proportional to 
S^-fj,j. To be specific, we use the Pauli-Sakurai 7 matrix convention in our calculation. 

We discretize our current operators by using the following relations for right and left derivatives 
in lattice [26] 
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Df,tp^{x) = — Uf,{x)i;^{x + a^)-Ul{x-af,)ij^{x-af,) 



2a 
1 

2a 



'4)^{x + a^) UUx) - ip^ix - a^) U^{x - a^) 



(9) 
(10) 



where a is the lattice spacing. For example, the two-index operator, On^, can be written as 



Of,y{x) 



A 



5a L 



i)''^\x)-i^U,{x)ij^f\x + a,) - ^^^\x)j^Ul{x-a,)^^f\x-a,) 



(/), 



+ V^^^^ (x - a,) 7/. U,{x - a,) V^^^^ (x) - ^^^^ (x + a,) 7^ Ul{x) (x) 
+ ^^^^ (x) 7, U^{x) (x + a^) - ^^^^ (x) 7, Ul{x - a^) V-^^^ (x - a^) 
+ ^^^^ (x - a^) 7, U^{x - a^) (x) - ^^^^ (x + a^) 7, ^7^ (x) V'^^^ (x)' 



(11) 
4. Similar 



where \ = —i for /U = 4, z/ = 1, 2, 3; A = +1 for = = 1, 2, 3 and A = —1 for /i = i 
expressions can be obtained for the three-index operators. 

Since in lattice QCD the continuous space-time space is described on a four-dimensional cubic 
lattice, the 0(4) group in the continuum reduces to the hyper-cubic group H{A) [27, 28]. This im- 
plies that operators belonging to irreducible representations of 0(4) may transform in a reducible 
way under H{A). This will allow them to mix with lower dimensional operators under renormaliza- 
tion. In order to avoid such mixing, it is suggested to adopt the following combination of operators 
which have minimal mixing (or no mixing) [29, 30, 31, 32]. For two-index operators (for (x)), we 



choose (i = 1,2,3) and O 
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O 



44 



-(Oil + O22 + O33) which does not suffer from any 



-{Oijj+OAkk) 



mixing [29, 31]. The best choice for three-index operator (for (x^)) is O/^u = O^u 

where i,j,k = 1, 2, 3( i ^ j 7^ A;), which still suffers from some mixing [29, 32]. The matrix elements 
for these operators are 



{P\Oi\P] 

{p\oL\p] 
{p\oL\p] 



2m 



{x)f^f-EP„ 



2m 



(x) 



— {x^) -Epl 
2m ■' ■' 



(12) 
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where the 2m factor, with m being the nucleon mass, is due to the normahzation of the spinors 
with u{p, s)u{p, s') = 6ss'- 



B. Two- Point and Three- Point Correlation Functions 



The proton two-point function we use (with the color indices suppressed) is 

G%%it,P} = Y.e-^P-^^--^^Ho\T[x''ix,t)x^ixo,to)m, (13) 

X 

where t is the nucleon sink time, and p is the momentum of the nucleon. The interpolating 
fields [34, 35, 36, 37, 38] we use are 

X^ix) = 6„fe,V'I(")'^(x)(C75)a/3<^'(x)^(")^(x), (14) 

Xyix) = -edefi^^:;}\x)^f\x){^^C),„i^^^^~^\x), (15) 

where u and d stand for up and down quarks, respectively. C = 7274, is the charge conjugation 
operator with the Pauli-Sakurai 7 matrices. The letters, a, 6, • • • , stand for the color indices. The 
Greek letters, a, /3, are the spin indices. 

Since we are interested only in nucleon with = , we use the projection operator T = 

1 TYl 

— (1 H 9374) [33] to eliminate the contamination from negative parity 5ii state. Here m~ and 

2 Ep 

Ep~ are the mass and energy of the state. After applying the projection operator, we get the 
two-point function as 

+ E E e-^^ (16) 

9=+- nW=l 

where the superscript 9 = +,— represents positive (negative) parity states and 

fin+'-,p) = iVr"^^(0|x"(xo)|n+'-,^?,s)(n+'-,p,s|x''(xo)|0), (17) 

s 

being the number of lattice points. 

As a result, the projected two-point function with momentum p at large time separation, i.e. 
t ^ to, will filter out the excited states, leaving only the nucleon state remaining asymptotically 

(t - to) > 1 



Tr [TGNN{t,p)] 



^-E°+(t-to)^ (18) 



where and Ep'^ are the nucleon mass and energy, respectively, k is the hopping parameter in 
the Wilson fermion action. From now on, we will drop the superscript +. 

The three-point function for any general operator O (with color indices suppressed) is defined 
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as 

X2,Xl 

where t = t2 is the nucleon sink time, t = ti is the current insertion time, t = to is the nucleon 
source time, and pi and pf are the initial and final momenta of the nucleon, respectively. For 
forward matrix element, pf = pi = p. In this case, 

G%%^(,t2,h,p) = J2 e-*^^-("°)(0|T(x"(f2,t2)0(xi,ti)x''(fo,io))|0). (20) 

X2,Xl 




(a) (b) 

FIG. 1: Quark line diagrams of the three-point function in the Euclidean path integral formalism, (a) 
Connected insertion and (b) disconnected insertion. 



The quark line connected part of the three-point function in the path integral is represented 
diagrammatically in Fig. 1(a). It needs to be stressed that it is not a Feynman diagram in pertur- 
bation theory. We see that the nucleon interpolating quark fields contract with the quark fields of 
the current so that the quark lines fiow continuously from t = to at the nucleon source to t = t2 at 
the nucleon sink. This is termed the connected insertion (C.I.). The quark line disconnected part is 
represented diagrammatically in Fig. 1(b). In this case, we see that the quark fields in the current 
self contract to form a loop, which is disconnected from the nucleon interpolating quark fields with 
regard to their quark lines. This is termed the disconnected insertion (D.I.). Since the quarks are 
propagating in the gauge background, if we were to consider a similar situation in terms of Feynman 
diagrams in the perturbative approach, it would involve gluon lines between the quark loop and 
the nucleon propagator so that the corresponding Feynman diagrams are connected in this sense. 
Indeed the corresponding disconnected Feynman diagrams are subtracted as the uncorrelated part 
in the definition of the disconnected insertion (D.I.) i.e. 

Gnon{B.I.) = ixOx) - {0){xx), (21) 
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where x is the nucleoli interpolation field. In the literature, this disconnected insertion is sometimes 
referred to as the "disconnected diagram" which can cause some confusion. 

The computation of the C.I. is relatively straight forward. We shall use the sequential source 
technique [30, 42, 43, 44] to calculate it. This fixes the source point to and the sink time slice t2- 
However, the computation of the D.I. poses a major numerical challenge. The D.I. contains not 
only the usual propagators from the source, xq, to any point, x, but also the propagators from any 
insertion position (xi) to any other lattice points. This amounts to inverting the fermion matrix at 
each point of the lattice to construct the all-to-all propagators. This entails inversion of a million 
by million (~ 16'^ x 24 x 3 x 4) sparse matrix on a 16^ x 24 lattice (3 and 4 are number of color and 
spin indices) for each gauge configuration. This is unattainable even by using the computing powers 
of today's supercomputers. Instead, we shall calculate with the stochastic method. Specifically, we 
adopt the complex Z2 noise for the estimation with unbiased subtraction. The detailed description 
of the method and the usefulness of discrete symmetries will be presented in Sec. Ill A. 

C. Ratios of Correlation functions 

In order to extract (x) and (x'^), we take the suitable ratios of the three-point to two-point 
correlation functions. Since, for the case of C.I., the nucleon sink time is fixed, and for the case of 
D.I., the sink time is not fixed, the procedures for extracting the matrix elements, after taking the 
ratios, are different for the two cases. 

1. Connected Insertions 

In view of the fact that the matrix elements of On, between the equal momentum nucleon states, 
are proportional to pi, and that of O/^a are proportional to pf, and that the momentum projection 
is folded in the sequential source at the sink time t2 in the connected insertion calculation, we 
have chosen only one momentum for the nucleon to be pi = 2tx jLa (the lowest available non-zero 
momentum) along the x direction in order to reduce the computational cost. 

After inserting complete sets of physical states between the interpolation fields and the operators 
O, we arrive at the asymptotic relations for the following ratios of three- to two-point functions at 
(ii - to) > 1 and (t2 - ti) > 1: 




(22) 



In practice, one takes a plateau in the insertion time, ti, to define the asymptotic region. 
Similarly, we get for the other two operators 



2kS^ 



Tr rG^(p^^^(t2,ti,pi) 



(ti - to) > l,(t2 -ti) > 1 




(23) 



T:[rG'iviv(t2,Pi)] 



and 



2k Tr rG^^^^^^(t2,ti,pi) {t 



to) > l,(t2 -ti) > 1 




(24) 



p1 Tr{YG^^{t2,pi)\ 
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2. Disconnected Insertions 



As seen from Fig. 1(b), the calculation of the valence quark propagators in the disconnected 
insertion is separate from the loop calculation in each configuration; this means that the momentum 
in the three-point function can be chosen in the nucleon two-point functions independent of the 
expensive loop calculation. For the present calculation, we choose pt = ±2-11 /La for the three-point 
function for the cases of and O/j^n, and zero momentum for O44. 

In contrast to C.I., the sink time need not be fixed in D.I. We can sum over the insertion time 
to gain more statistics [56, 57, 58, 60]. There are various methods by which the summation is 
performed [47, 60, 66, 72]. All these methods, except in [60], considered summation over insertion 
time up to or beyond sink time. From Appendix A, we see that the matrix elements for the 
twist-two operators are analytically zero if ti < Iq and ti > t2- (It may not be zero for other 
operators. Even then it will not contribute to the physical quantity intended to be measured). So 
the summation outside the nucleon source and sink times will contribute to unnecessary noise and 
possible contribution from higher states which are unrelated to the target matrix element. Also, at 
the source and sink time it can have contributions from the contact term. In view of the above, we 
shall take the sum to be from (source time -|-1) to (sink time —1) [60]. According to the derivations 
in Appendix A 2 (see Eq. (AlO)), under the condition t2 — to ^ 1, we get for the operator On 



t2-l 

E 

tl=to + l 



2k Tr [GN04iN{t2,ti,P'i)] 
Pi Tj:[T GNN{t2,Pi)] 



{t2 - to) > 1 



{x)ii t2 + const. 



(25) 



Similarly, for the other two operators, we get 

'^7VO44^(*2,tl,0) (t2-io)>l 



2k ^ 



(x)44 12 + const.. 



(26) 



and 



E 



'^7V04,,A^(*2,tl,K) (t2-to)»l 



t,=t,+iP'^ Tr[rG^^(t2,K)] 



{x'^)Aiit2 + const. 



(27) 



III. DISCONNECTED INSERTION CALCULATION 



The D.I. calculation is the most numerically intensive part. We shall discuss the various aspects 
of the calculation in more detail. 



A. Discrete Symmetries and Transformations 



Since the D.I. calculations are performed by using a stochastic noise estimator, the signals for the 
current loop are always noisy. We can reduce the errors by making good use of some of the discrete 
symmetries, specifically parity, 75 hermiticity, and charge-75 hermiticity {CH transformation) [45, 
46, 47, 61]. 
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By applying these symmetries and transformations, one can then work out the effects of gauge 
averaging in order to find out the correct part (e.g. even or odd parity, real or imaginary part, 
etc.) of the two-point and three-point correlation functions, current operators etc., and discard the 
irrelevant part. 

1. Two-point Functions 

Since in the case of D.I. the three-point functions are constructed by multiplying the nucleon 
propagator with the current loop in each gauge configuration (Fig. 1(b)), we have the advantage 
of taking into account nucleons propagators with equal and opposite momenta in order to increase 
statistics. While doing so, we have to consider the fact that such combinations have appropriate 
parity and have appropriate real or imaginary part w.r.t. the loop. In Table I, we show the effect 
of parity and CH transformations on various such combinations. We denote Tr \r GNNit^P] U)] to 
be the nucleon propagator on each gauge-configuration with momentum p. 



Nucleon Propagators 


CH 

Transfer mat ions 


Parity 


TT[TGNN{t,p-.U)] 


TT[TGNN{t,-p;U*)Y 


T,[TGNN{t,-p\UP)] 


Tr [TGNN{t,p;U)] 
-Tr [TGNN{t,-p;U)] 


-{Tr [TGNN{t,P\U*)] 
-TV [rGiviv(t,-p; [/*)]}* 


Odd 


Tr [TGNN{t,p-,U)] 
+Tr [TGNN{t,-p;U)] 


{tV [TGNN{t,P]U*)] 
+TY [TGNN{t,-p;U*)]Y 


Even 



TABLE I: Table showing the outcome of the parity and CH and parity transformations on the combinations 
of nucleon propagators with equal and opposite momenta. U'^ is the parity transformed gauge link. 



2. Current Loop 

The outcome of the parity, 75 hermiticity, and CH transformations for the operators On, ©44, 
and O/iii on each gauge configuration are shown in Table II. The notation "Im" includes the factor 
i of the imaginary part of each operator. 

3. Correlations between the Nucleon Propagator and the Loop 

When determining the correlation between the appropriate parts of the nucleon propagator and 
the loop, we have to consider the three-point function as a whole. As an example, if we consider 
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Loops 


Parity 


75 

Hermiticity 


CH 

Transfor mat ions 


75 Hermiticity 
CH Transformations 
Combined 


Oii 


Odd 


Imaginary 




Im [L]^ = Im [L]^, 


-^(011 + 022 + 033) 


Even 


Real 




Re[L]^ = Re[L]f,. 


O4M — \ {Oijj + Oikk) 


Even 


Imaginary 




Im [V\^ = -Im [V\^, 



TABLE II: Table showing the outcome of the parity transformations, 75 Hermiticity and CH Transforma- 
tions on current loops. L stands for current loop and "Im" includes the factor i of the imaginary part. 



two-point functions with both the momenta ±pj for the operator O^a, then from Eq. (25), we get 

y ^ [GN0Mt,M,V^)] -^^ [GNO,Mt2M.-p.)] ^ + (28) 

So, for which is odd in parity, one can explicitly consider the odd combination of the nucleon 
propagators in Table I, and similarly, the even combinations for ©44 and Om, to reduce noise. 
According to the CH theorem [45, 46, 47], the path integral for {O) in QCD is either real or 
imaginary (except in the case with chemical potential). It can be shown from Tables I and II that 
if we combine the CH transformation of the nucleon propagators (with the appropriate combination 
or parity) and the loop, the three-point functions transform to the positive complex conjugate of 
themselves on the link U for all the three D.I. cases considered here. This, according to the CH 
theorem, means that the three-point functions are real. Since, by using 75 hermiticity we see that 
the loops are either real or imaginary (Table II); one needs to multiply them with only the real or 
imaginary part of the nucleon propagators to make the three-point functions real. For example, 
the L.H.S. of Eq. (28) can be written as 

^ [Tr [GNO^^Nih^ti^Pi)] -Tr [GNOuN{t2,ti,-pi)]] 



— E 



^ [TV [rGNN{t,pi;W)] - TV [V GNN{t, -Pi;W)]] x Im(Loop)^i, 



+ ^ [Tr [TGNN{t,pf,U*^)] - Tr [T GNN{t, -pi;U*^)]] x Im(Loop)c;,, 

= Im {IV [rGNN{t,Pi;W)] - Tr [V GNN{t, -p^■,W)]} x Im(Loop)t;.. (29) 

^ j 

where Ng is the number of gauge configurations. As illustrated above, one can exclude the real 
part of the nucleon propagator and the real part of the loop for O^i operator which contributes to 
noise with finite number of noise vectors. The utilization of 75 hermiticity and CH theorem has 
been shown to reduce noise effectively in the loop calculation of the quark angular momentum [47] . 
Similar procedures can be applied to the other two operators. We show in Table III the relevant 
parts of the nucleon propagator(s) and the corresponding loop in each row which are to be correlated 
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on each gauge configuration. 



Nucleon Propagators 


Loops 


Im[Tr [T GNNit,p;U)] - Tt [T GNNit,-p;U)]] 


Im[04i] 


Re [TV [rGNN{t,0;U)]] 


Re [O44 - ^ (On + O22 + O33)] 


Im [i [Tr [^GNN{t,P^;U)]+TT [F GNN{t, -PuU)]]] 


Im [Oi^^ - \ (O^jj + 04fcfe)] 



TABLE III: The relevant parts considered for DT. calculations for the nucleon propagator and the corre- 
sponding loops based on parity, 75 hermiticity, and CH theorem. 



B. Complex Zi Noise and Unbiased Subtraction Method 

As we mentioned in Sec. II B, it is a numerical challenge to evaluate the quark loop. We shall 
adopt the complex Z2 noise with unbiased subtraction to calculate it. 

The basic idea of the complex Z2 noise method [51] is to construct L noise vectors, r?^, 77^ 77^ 

(each of dimension x 1), where rf = {ry-J, r/2, •.••t?;^}, in order to stochastically estimate the 

1 

inversion of an A^ x A^ matrix. Each element r/^ takes one of the four values, — p{ibl it i}, chosen 

v2 

independently with equal probability. They have the properties of a white noise 

{rii) = hm 7 V ??r = 0. ivjvj) = 1™ y V ??J"r/" = 6ij, ri\r]i = 1. (30) 

L—*oo L — ' L^cx) L ^ — ' 

n=l n=l 

Then, the expectation value of the matrix element M~j^ is obtained by solving for Xi in the matrix 
equations MX = rj, so that 

E[Mri] = (r/Jx.) =Y^M-Hri^^ri,) = M7^\ (31) 
k 

It has been shown [51, 52, 53] that the variance corresponding to the estimator is given by 

(32) 

which is minimal, since it does not involve the positive contribution from the diagonal matrix 
elements as do in other type of noises, such as the Gaussian noise. It is in this sense, Z2 noise 
is considered optimal [51, 52, 53]. Similarly, one can show that Z{N) and U{\) noises are also 
optimal. 

It has been further shown that the off-diagonal matrix element contributions to the variance 
in Eq. (32) can be reduced by subtracting a judiciously chosen set of traceless N x N matrices 

N 

Qip) [47^ 53^ 61]^ which satisfy Qnji = 0; P = 1; " " " i -P- Then the expectation value is unchanged 

n=l 



2 



-1 
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when M ^ is substituted with M ^ — Xp Q^^^ (Ap is a constant) 



E[(r?t(M-i-^ApQ(P))r/)] = Tv M-\ 
p=i 

while the variance becomes 



(33) 



4^ = Var[(j7t (M 



-1 



P N 



M~ 



iv) 

m,n 



(34) 



So, with a judicious choice of traceless matrix, the variance may be reduced when the ofF-diagonal 
matrix elements of are correlated with those of M~^. This subtraction is unbiased, because 
it doesn't change the expectation value of the trace. The natural choice for the set of traceless 
matrices is the hopping parameter expansion of the inverse of the Wilson fermion matrix, M [53] , 
given by 



M"^ = I + kD + k'^D'^ + K^D^ + 



where 



Dx,y = E (l-7/.)a/3C/';''(x)<5^,y-a^ + {1 + 7^,)aP Ul"^ {x - a^) 6^^y+a, 

,1=1 L 

which is off-diagonal in space-time. 

In evaluating the quark loop, we shall consider the following subtraction 

( 7?t Q' [M-i -{I + kD + k^D^ + K-^D^ + ...)] r] ) + traces, 



(35) 



(36) 



(37) 



where Q' represent operators between the fermion fields in 0^,0^4,, and Om. Traces may arise 
while combining the operators Q' with the hopping expansion in Eq. (35), leading to non-zero traces 
at certain power D. In this case, one needs to evaluate these traces and add them back to cancel 
out those which are subtracted out in the noise estimation. The details of the trace calculations for 
these operators are given in Appendix B. Since the calculations of traces are cumbersome for large 
loops, we will restrict ourselves only up to k^L)^ for two-index operators and k'^D^ for three-index 
operators. We shall now discuss these trace terms here. 



Traces for Two-Index Operators O^i and On 



Since Q' for On and are point-split, the trace for the first subtraction term, i.e. Tr((5' I), is 
zero. So we can use it for subtraction. However, from Appendix B, we see that the noise estimate 
is real for O^i while the exact loop Tic{Q' M~^) is imaginary (see Sec. IIIA). Similarly, the noise 
estimate is imaginary for O^^ while the exact loop Tv{Q' M~^) is real. Therefore, there is no use to 
include this subtraction term. The subtraction with k^D^ and k'^D'^ terms in the hopping expansion 
are also traceless, since the multiplication of Q' to these terms does not lead to a plaquette which 
is the lowest order in k that contributes to a trace. On the other hand, the kD and k^D"^ terms 
can lead to a plaquette, thus can have non-zero traces. Since the operators, O^i and 0^fj_, lead to 
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different traces, we sliall consider them separately. 

In the case of O^i, the contribution from the hopping expansion terms kD, k^D'^, k'^D^ and 
k'^D'^ are traceless (Appendix B). One does not have to worry about the complication of having to 
add back the trace contributions. 







1 i 1 

X + + a^/ X — + aj^ 


4 




3rd term 


1st term 


^^ 


2nd term 




4th term 


X 

2nd term 




X 

1st term 










a; + 04 - 



(a) 



(b) 



FIG. 2: Plaquette terms (a) for the operator when k D term is considered and (b) for Om when n D 
term is considered. 

But, operator has non-zero traces when combined with kD and k^D^ (see Appendix B for 
a full derivation). The trace for the kD term is — OGkVs, where V3 is the three- volume of the lattice 
used, and A = -|-1 for = = 1, 2, 3 and A = —1 for jjL = v = A. And, for tz^D^ [see Fig. 2(a)], it is 

+ U^{xi)ul^,{xi +a^- a^')C/^(xi - a^>)U^>{xi - a^>) 
+ U^{xi -af, + a^/)C/^,(xi)C/^(xi - a^)C/^'(xi - a^) 

+ U^i{xi - a^/)C/^(xi - a^,)Ul,{xi - - a^/)C/^(xi - a^, - a^/)}™"" . (38) 
These traces need to be added as shown in Eq. (37). 



2. Traces for Three-Index Operators 



For the three-index operator, the noise estimator is real for the first term, I. But from 75 
hermiticity we see that the loop for 044^ is imaginary. Thus, it is not useful to consider this term 
in subtraction. Up to k^D^, the only trace contribution comes from the k^D^ term, which is (see 
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Appendix B for a full derivation) 

2^16k^ ^ ^lm{ul{xi - ai)Ui{xi - ai)U4{xi)Uj {xi + 04 - Oj)}^ 



+ ^lm{U4{xi - 04 + ai)U- {xi)ul{xi - a4)Ui{xi - 04)}'' 

m 

The graphical representation for this trace is illustrated in Fig. 2(b). 



(39) 



3. Numerical Test 



It is necessary to check that the analytical calculations for trace for each subtraction term is 
correct. In order to do that we have to calculate the subtraction terms exactly. We did it on a 
smaller (4^) lattice as described below. 

From the matrix equation MX = rj, we have 



= (I + kD + K^D^ + tC-D^ + K^D^ + ■ ■■)r]r{^ . 



(40) 



From Eq. (40) we see that we can exactly calculate M ^ i.e. each subtraction term, only if r]ri^ is 
a unit matrix. By choosing the following orthogonal set of vectors {rji\ 



/1\ 







1 







1 



VJ VJ VJ 






VI/ 



(41) 



1 

we obtain rji]^ = y ??i t]} = 1- Using these vectors, we compute I, kD, k'^D'^, ■ ■ ■ for each oper- 

i 

ator. We have found that our numerical results matched with those of analytical expressions e.g. 
for O/^i operator, and kD, k^D'^, k^D^ and k'^D^ terms have no traces. 



IV. RENORMALIZATION 



The physical matrix elements, which are determined from a linear extrapolation to the chiral 
limit, are extracted in lattice units from the Monte Carlo calculation. In order to relate to the 
experimental values, they have to be expressed in physical units and renormalized at a certain 
scale. The renormalized operators, at a finite energy scale //, are related to the bare lattice 
operators through the renormalization constant 

Oifi) = Zo{a^,,g{a))0{a), (42) 
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where g is the bare couphng constant, which is equal to one in our case. The renormahzation con- 
stants, Z's, are computed by using perturbation theory. We will use the values of Z factors, which 
are computed in [62] by using tadpole improved perturbation theory [64]. Since the experimental 
results are often renormalized in the MS scheme, we will also use the calculated Z factors matched 
to this scheme. In the quenched approximation, the renormahzation constants are 



Zo{an,g*) 



Up 



9 



*2 



Cf ( 70 ln(a/i) + B^^ + (no - 1) Stt^Zq ) + O {g*^) ) , (43) 



where 



0.155, 



(44) 







and 70 is the anomalous dimension of the operator, Bo is the finite part of Zo, and n/j is the 

number of covariant derivative(s) in the operator. We list 70 and Bq in Table IV for the three 
operators we consider. 



Operators 


10 


rMS 


O^i 


16 


1.279 


Oaa 


16 
3 


2.561 




25 
3 


-12.128 



TABLE IV: The values of 70 and _B^^ for all the three operators under consideration. 

Our inverse lattice spacing is determined to be 1/a = 1.74 GeV [47] by using nucleon mass. The 
values of B^^'s are given in Table IV. The value of uq = (| Tr t^pi^q)"' ~ ^-^^ obtained from 
Ref. [59]. By using all the relevant factors for a particular operator, we get the renormahzation 
factors at /U = 2 GeV (^^ = 4GeV2) scale for the following operators as 



04i : 




= 0.972, 


O44: 


^044 


= 0.953, 


Oah : 




= 1.116. 



(45) 



Also, in the tadpole-improved mean-field approach, there is a finite ma correction factor / for 
the Wilson fermion in the case of fermion bilinear operators [57, 59, 65] 



/ 



1 



m„a 



ln^-3 , 

K 



(46) 



with the critical Ate = 0.1568 [59]. The values of / for k = 0.154,0.155, and 0.1555 (which are used 
for our calculation) are 0.972, 0.948, and 0.936, respectively. 
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V. NUMERICAL PARAMETERS AND ERROR STUDIES 

We use 500 gauge configurations on a 16^ x 24 lattice generated with Wilson action at /3 = 6.0 
in the quenched approximation. They are produced by the pseudo-heatbath algorithm with 10,000 
sweeps between consecutive configurations. The values of the hopping parameter we have used 
are k = 0.154, 0.155 and 0.1555. The critical hopping parameter, = 0.1568 is obtained by a 
linear extrapolation to the zero pion mass [59]. Using the nucleon mass to set the lattice spacing 
at a = 0.11 fm, the corresponding pion masses are 650(3), 538(4), and 478(4) MeV, and the 
nucleon masses are 1291(9), 1159(11), and 1093(13) MeV, respectively. We used Dirichlet boundary 
condition in the present work. We should note that there is a large uncertainty in determining the 
lattice spacing in the quenched approximation, as much as ~ 20%. For example, using ro = 0.5 fm 
to set the scale, the lattice spacing would be a = 0.09 fm. Thus, using this scale, the dimensionful 
quantities will be shifted by ~ 20%. Since we are calculating the moments (x) and (x^) which 
are dimensionless, the results we report here will not depend on the scale of the lattice spacing 
except for the renormalization constant which has an negligible difference between these two lattice 
spacings. 

• For the connected insertion, we have chosen the number of independently generated gauge 
configurations to be 200 and for disconnected insertion it is 500. The maximum number of 
noise used is 400 for each gauge configuration. 

• The error analysis has been performed by using the jackknife procedure [67, 68]. Since the 
computations for all the time slices and quark masses has been performed using the same 
set of gauge configurations, the data we obtain are correlated in Euclidean time. The corre- 
lation among different quantities are taken into account by constructing the corresponding 
covariance matrices [69, 70, 71]. The error bars are obtained by using this method. In order 
to extract various physical quantities, we have used correlated least-x^ fits. 

A. Error Studies in DI 

Since noise estimate plays an essential role for the case of the disconnected insertion, we shall 
study the effect of number of noise and gauge configurations on the error for the signal we are 
extracting. We will present a few results of such studies before discussing the results for the moments 
of the quark distribution. 

The standard error due to the Z2 noise estimation of the loop averaged over the gauge config- 
urations is given by [54, 73, 74] 



the noise estimator, and Ng and Nn are number of gauge configurations and noise estimator re- 
spectively. 

To make sure that our results are generated correctly, we verify that their errors are in conformity 
with Eq. (47) and extract the standard deviations ag and (T„ from the existing data with and without 
unbiased subtractions.. 




(47) 



where ag 



is the variance of the ensemble of gauge configurations, cj^ is the variance coming from 
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FIG. 3: Errors of the noise estimation plotted against the number of configurations for different sets of noise 
vectors for the loop part of the current, O^i at Kg — 0.154 nd insertion time, ti = 14 (a) without subtraction 
and (b) with four subtraction terms. 
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FIG. 4: Errors of the noise estimation plotted against the number of noise vectors for different sets of 
configurations for the loop part of the current O/^i at Ks = 0.154 and insertion time, ti = 14 (a) without 
subtraction and (b) with four subtraction terms. 



In Figs. 3(a) and 3(b), we plot the errors of the noise estimation for the loop part of the current 
in Eq. (47) at the time slice 14 against the number of gauge configurations for 30, 100 and 
400 noise without unbiased subtraction and with four subtraction terms {kD, k^D'^, k^D^ and 
k'^D^) respectively. The errors on the errors are obtained by using the double jackknife method 
on the data. After fitting for and from Eq. (47), we see that all the curves can be well 
described by Eq. (47). Similarly, we plot the errors against the number of noises for 100, 250 and 
500 gauge configurations without subtraction and with four subtraction terms in Figs. 4(a) and 4(b) 
respectively. Again, we see that the curves fit Eq. (47) well. A similar conclusion can be drawn for 
the other two operators. The central values of the standard deviations, ag and an along with their 
errors are given in Table V for the cases with and without subtractions. These values show that 
the standard deviations for the gauge configuration and noise ensemble are not of the same order. 
The standard deviation for noise are, in fact, much higher than that for the gauge configuration. 
Another point to note is that the standard deviation an for the noise is reduced by almost a half 
with 4-term unbiased subtraction whereas ag remains the same. 



16 





'-'^gauge 


•^noise 


No Subtraction 


0.204 ± 0.063 


9.341 ± 0.301 


With 4 Subtractions 


0.205 ± 0.028 


5.500 ± 0.201 



TABLE V: Table for the values of standard deviations of gauge configurations and noise for the current Om 
without subtraction and with four subtraction terms. 



Methods 


Fitting Range 


No Subtraction 


4 Subtractions 


Summation from 5 to 20 


10 - 14 


0.011 ± 0.019 


0.030 ± 0.012 


Summation from source to sink time 


10 - 14 


0.007 ± 0.013 


0.029 ± 0.009 


R. Lewis et al. [72] 


11 - 14 


0.000 ± 0.012 


0.027 ± 0.009 


Summation from (source time + 1) 
to (sink time + 1) 


10 - 14 


-0.002 ± 0.014 


0.026 ± 0.010 


Current analysis 


10 - 14 


0.004 ± 0.012 


0.028 ± 0.008 



TABLE VI: Table for the vales of {x) (D.I.) at = Kg = 0.154 for the operator by using five different 
methods. 



B. Analysis for DI 

As mentioned in Sec. II C 2, we have studied five different methods of summation over insertion 
time by using the operator at = Ks = 0.154 shown in Fig. 5. In the first method [66], we 
have performed the summation of the current insertion starting and ending 4 time slices away from 
each of the boundary. In our case, it would be from 5 to 20 [Fig. 5(a)]. In the second method [47], 
the summation has been performed from the source to the sink time of the nucleon propagator 
[Fig. 5(b)]. The third method is described in [72] [Fig. 5(c)]. The fourth method is an additional 
study where the summation has been performed from (source time + 1) to (sink time + 1) of the 
nucleon propagator [Fig. 5(d)]. The fifth method used in [60] and described in Sec. II C 2 [Fig. 6(a)]. 
In this method, the summation has been performed from (source time + 1) to (sink time — 1) 
of the nucleon propagator. In the first, second, fourth, and fifth methods, the slopes (given by 
Eq. (25)) are fitted between the time slices 10 and 14 in order to extract the signal. And in the 
third method, a constant is fitted between 11 and 14. The values are provided in Table VI. We 
see that the all these methods are consistent with each other. For our present work, we adopt the 
fifth method. 

C. Multiple Sources for DI 

In the case of D.I., we have the liberty of choosing quark propagators at different source loca- 
tions on the lattice with much less overhead that the C.I. computation which involves generating 
additional quark propagators using the sequential source method. We correlate these nucleon 
propagators at different source locations with the already computed loop to increase statistics. We 
have used 1, 4, and 16 different sources. This results in significant reduction of error bars which is 
presented in Sec. VI. 
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FIG. 5: The ratio (D.I.) of the three-point to two-point functions at k^, = Ks — 0.154 for the O^i operator 
is plotted against the nucleon sink time (<2) by using four different methods: (a) summation of insertion 
time from 5 to 20, (b) summation of insertion time from source to sink time, (c) method used in R. Lewis 
et al. [72], and (d) summation of insertion time from (source time -I- 1) to (sink time + 1). 



VI. RESULTS 



In this section, we will present results for the first and second moments of nucleon's parton 
distribution function for both the disconnected insertion and the connected insertion. 



A. Disconnected Insertions 



The disconnected insertion [Fig. 1(b)] entails the correlation between the quark loops with up, 
down, and strange quark currents and the nucleon two-point propagators as discussed in Sec. II B. 
We consider the up and down to have the same mass so that their moments are the same. 



1. First Moments 



First, we will discuss the results for the first moments from the 04i(i = 1,2,3) operator for 
the case where the valence quarks in the nucleon propagators and the sea quark in the loop are 
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the same. This is the case, when extrapolated to the physical u/d quark mass would give the 
disconnected insertion result for the u and d quarks. We average over the three spatial directions 
on each configuration. From now on, Oin stands for the average over 1, 2 and 3 directions. 
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FIG. 6: The ratio of the three-point to two-point functions (summed over insertion time) for the 
operator, for the case with equal valence and sea quark masses, is plotted against the nucleon sink time (^2) 
at (a) Ky — Ks — 0.154, (b) — Ks = 0.155, and (c) = Kg — 0.1555. (d) is a linear extrapolation to the 
chiral limit for the first moment, {x)u,d, of the up (down) quark which is plotted against niqa. 



For the disconnected insertion, we shall define two k's for the quark mass: for valence 
quarks, and Ksea for sea quarks. For the strange quark currents we have fixed Ksea = 0.154, which 
is close to the strange quark mass as determined from the (/> meson mass and takes the values 
of 0.154, 0.155, and 0.1555. We consider the cases with equal valence and sea quark masses, 
i.e. Ksea = = 0.154,0.155, and 0.1555 in order to extrapolate to the chiral limit to obtain 
{x)u+u = {x)d+d = {x)u4 (D.I.). 

In Figs. 6(a), 6(b), and 6(c) we plot the ratios in Eq. (25) against the nucleon sink time t2- The 
insertion time is summed from [source time -|- 1] (i.e. to + 1) to [sink time — 1] (i.e. t2 — 1). In 
these figures, we see that, after the unbiased subtraction of four terms, there is a clear straight line 
behavior starting from the time slice 10. On the other hand, without the unbiased subtraction, 
there is not a clear signal. In fact, they are consistent with zero slopes. Also, the plots show that 
the error bars get reduced after subtraction by a factor of ~ 1.5. To extract the values of {x)u,d at 
each Ky, we have performed a correlated fit of the slope between the time slices 10 and 14. It gives 
us the value of the (x) at the corresponding Ky. The values of {x)u,d (D-I-)) along with their errors. 
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FIG. 7: The ratio of the three-point to two-point functions for the O^i operator is plotted against the nucleon 
sink time (^2) at Ky — 0.154 and Ks — 0.154 for 1, 4, and 16 sources after 4 subtractions. 
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FIG. 8: Linear extrapolation to the chiral limit for the first moment with 1, 4, and 16 nucleon sources (a) 
for {x)u,d (D.I.) and (b) for (x), (D.I.). 



are listed in Table VII. 

In chiral perturbation theory, the first moment of (x) has a leading non-analytic behavior 
oc m^ln(m^//i^) and leading analytic behavior oc [75, 76, 77, 78]. Since our pion masses are 
relatively heavy, we do not expect to be in the region where the non-analytic behavior is important. 
Furthermore, our present calculation is based on the quenched approximation. In view of this and 
other systematic errors that we have not taken into account, such as the large volume limit and 
continuum limit, we shall take the conservative linear extrapolation of {x)u,d to the chiral limit 
with the form A + Bruqa [Fig. 6(d)]. This linear extrapolation will inevitably introduce systematic 
errors. This issue will be dealt with when the configurations with lighter quark masses are available. 
We will certainly include the non-analytic behavior when the results from the PAC-CS 2-|-l-flavor 
dynamical clover fermion configurations are available with lighter sea quark masses [81]. Before 
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FIG. 9: The ratio of the three-point to two-point functions (summed over insertion time) for the O^i 
operator, for the strange quark which is iixed a.t k = 0.154, is plotted against the nucleon sink time (^2) at 
(a) — 0.154 and Kg = 0.154, (b) = 0.155 and Kg = 0.154, and (c) Ki, = 0.1555 and Kg = 0.154. (d) is a 
linear extrapolation of the valence quarks to the chiral limit for the first moment, (x)s, of the strange quark 
which is plotted against (m -I- ms)a. 



extrapolation, we have converted the values of {x)u^d (D.I.) to those of tadpole improved values 
by using the factors in Eq. (46). We did the similar analysis for the cases with 4 and 16 nucleon 
sources. Fig. 7 shows the ratios between three-point to two-point functions for all the sources with 
4 subtractions, and Fig. 8(a) shows the linear extrapolation to the chiral limit. In Table VII, we 
list the renormalized and linearly extrapolated values (to the chiral limit) of {x)u^d along with their 
errors. As stated in Eq. (45), the renormalization factor for this operator is 0.972. We find that the 
values for {x)u,d are about 3a away from zero for one and four sources and more than 5a for 16 
sources. We consider this a solid affirmation that we have been able to calculate the D.I. of {x)u,d 
via the noise method. 

Next, we consider the strange quark loop. This time we have fixed Ksea = 0.154 which corre- 
sponds to the strange quark mass. In Figs. 9(a), 9(b) and 9(c), we have plotted the ratio in Eq. (25) 
against the nucleon sink time by using the valence quark masses at = 0.154,0.155, and 0.1555, 
respectively. A similar procedure has been followed as for the up (down) currents to obtain the 
slopes, which gives the values of {x)s- Again, from these figures, it is clear that if we do not use un- 
biased subtraction, we will not see a signal. The error bars get reduced after subtraction. The linear 
extrapolation to the chiral limit is performed with the form A + B{rn + ms) where [m + ms) oc rn'j^, 
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1 source 
(No sub) 


1 source 
(4 sub) 


4 sources 
(4 sub) 


16 sources 
(4 sub) 


{x)u,d 

(D.i.j 


Ky = Ksea — 0.154 


0.004 ± 0.012 


0.028 ± 0.008 


0.013 ± 0.005 


0.016 ± 0.003 


= Ksea — 0.155 


0.009 ± 0.017 


0.037 ± 0.012 


0.018 ± 0.007 


0.022 ± 0.004 


= Ksea — 0.1555 


0.011 ± 0.022 


0.044 ± 0.014 


0.021 ± 0.008 


0.025 ± 0.005 


Linear Extrapolation 


0.033 ± 0.058 


0.056 ± 0.019 


0.028 ± 0.011 


0.032 ± 0.006 


{x)s+s 

(D.I.) 


= 0.154, Ksea = 0.154 


0.004 ± 0.012 


0.028 ± 0.008 


0.013 ± 0.005 


0.016 ± 0.003 


Kjj = 0.155, Ksea = 0.154 


0.005 ± 0.017 


0.034 ± 0.011 


0.016 ± 0.006 


0.020 ± 0.004 


= 0.1555, Ksea = 0.154 


0.004 ± 0.021 


0.038 ± 0.014 


0.017 ± 0.008 


0.023 ± 0.004 


Linear Extrapolation 


0.005 ± 0.029 


0.046 ± 0.018 


0.021 ± 0.010 


0.027 ± 0.006 


{x)s+s 




0.69 ± 0.64 


0.85 ± 0.13 


0.95 ± 0.18 


0.88 ± 0.07 





TABLE VII: (x) (D.I.) for up (down) and strange quarks at various k's and the linearly extrapolated results 
to the chiral limit with different number of nucleon sources for O^i operator. 



with fn being the average of the up and down quark masses and being the strange mass. In 
Table VII, we have listed values of {x)s+s (D-I-) for different Kt,'s and fixed Ksea- Also listed are 
the linearly extrapolated values of {x)s+s (D-I-) to the chiral limit. Again, we find that the value 
for {x)s+s is 4.5(7 away from zero for 16 sources. We have also listed the ratio of {x)s+s (D-I-) 
to {x)u+u (D-I-) in Table VII. We see that the ratio is close to 1 as expected [21]. We should 
point out the caveat that these results are based on the linear extrapolation to the chiral limit 
which are subjected to systematic corrections as mentioned above. The ratio is expected to be less 
susceptible to the systematic errors except the chiral extrapolation due to the non-analytic terms. 
In view of the fact that the lattice calculations of both the {x)u+d and {x)u-d for the connected 
insertions at quark masses between the strange and the physical u/d mass are all larger than the 
respective experimental results (see Table XII in Sec. VI B), we expect the ratio of {x)s+s (D-I-) to 
{x)u+u (D.I.) for the physical u/d mass to be larger than that obtained with a linear extrapolation 
to the chiral limit. 

We have performed the same analysis for the — | (On + O22 + C33) operator as in the case 
of the Oiii operator. The only difference is that we have used the nucleon with zero momentum in 
this case. But, due to the subtraction of the spatial trace terms in the current, it has significant 
numerical cancelations resulting in large statistical errors [29, 31]. The error bar of the slope is ~ 
3-7 times larger than that of Oin operator. Due to this large error, the signal is only ~ la to 
1.5cj away from zero shown in Table VIII. This happens for both {x)s+s and {x)u^d (D.I.) for this 
operator. 

2. Second Moments 

We will now present the results for the second moments obtained by using the current Oim — 
^ {Oijj + O^kk) {i ¥^ j ¥^ k; i,j,k = 1,2,3) for up (down) and strange quarks (see Figs. 10, 11, 
12 and 13. We average over the results from the three operators: O411 — | (©422 + O433), O422 — 
\ (C433 + O411) and 0433 - \ (0411 + 0422)- So O/^u - \ (O^jj + O^kk) will mean average over 1, 
2 and 3 directions from now on. In addition, we have used three subtraction terms (k-D, k^D^, 
and K?D^). For this operator, we have performed the same analysis as for the operators for the 
first moments, with the fitting done from t2 = 11 to 14. The values are presented in Table IX. 
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1 source 


1 source 


4 sources 


16 sources 






(No sub) 


(4 sub) 


(4 sub) 


(4 sub) 




= K,sea — 0.154 


0.045 ± 0.044 


0.033 ± 0.040 


0.025 ± 0.026 


0.038 ± 0.020 


{x)u,d 


— l^sea — 0.155 


0.075 ± 0.056 


0.059 ± 0.053 


0.032 ± 0.035 


0.046 ± 0.027 


(D.I.) 


= Ksea — 0.1555 


0.095 ± 0.065 


0.073 ± 0.060 


0.034 ± 0.042 


0.049 ± 0.033 


Linear Extrapolation 


0.130 ± 0.081 


0.095 ± 0.071 


0.037 ± 0.050 


0.058 ± 0.043 




K„ = 0.154, Ksea = 0.154 


0.045 ± 0.044 


0.033 ± 0.040 


0.025 ± 0.026 


0.038 ± 0.020 


{x)s+s 


K„ = 0.155, Ksea = 0.154 


0.067 ± 0.056 


0.047 ± 0.050 


0.026 ± 0.034 


0.041 ± 0.027 


(D.I.) 


= 0.1555, Ksea = 0.154 


0.087 ± 0.065 


0.062 ± 0.057 


0.026 ± 0.039 


0.041 ± 0.033 


Linear Extrapolation 


0.114 ± 0.080 


0.077 ± 0.068 


0.026 ± 0.048 


0.043 ± 0.042 



TABLE VIII: Table for the values of (x) (D.I.) for up (down) and strange quarks at various kappa values 
and after linear extrapolation to the chiral limit with different number of nucleon sources for the ©44 — 
i (On + O22 + C33) operator. 



No Subtraction ' 
3 Subtractions i 



500 Configurations and 400 Noise| 



Summation of Insertion Time: (source time+1) to (sink time-1) 



Slope after 3 subtractions (One Source)= 0.0004 ± Q.0035 | 



10 12 
sink time (t2) 



(a) 



I No Subtraction 
I 3 Subtractions ' 



500 Configurations and 400 Noisel 

Summation of Insertion Time: (source time+1) to (sink time-l)| 
Slope after 3 subtractions (One Source)= 0.0104 ± 0.0095 | 



10 12 
sink time (12, 



(c) 







No Subtraction •■ 
3 Subtractions i 



500 Configurations and 400 Noise] 

Summation of Insertion Time: (source time+1) to (sink time-l)| 



0.15 Slope after 3 subtractions (One Source)= Q.0045 ± 0.0065 | 



10 12 
sink time (t2) 



(b) 



0.035 I 

0.03 
0.025 

0.02 - 
0.015 

0.01 

0.005 - 

~ - 

-0.005 

-0.01 



Chiral extrapolation for <x j| 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 

a 



(d) 



FIG. 10: The ratio of the three-point to two-point functions (summed over insertion time) of the 04^ — 
^ {04jj + Oikk) operator, for up (down), is plotted against the nucleon sink time ^2 at (a) k„ = Ks = 0.154, 
(b) Kv = Ks = 0.155, and (c) — Ks — 0.1555. (d) is the linear extrapolation to the chiral limit for the 
second moment, {x'^)u,d (D-I-) plotted against niqa. 
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0.04 
0.02 


-0.02 
-0.04 
-0.06 
-0.08 
-0.1 



k,=0.154, Kjja=0.154 



No Subtraction ' 
3 Subtractions i 



500 Configurations and 400 Noise| 



Summation of Insertion Time: (source time+1) to (sink time-1) 



Slope after 3 subtractions (One Source)= Q.0Q04 ± 0.0035 | 



10 12 
sink time (t2) 



(a) 



0.15 
0.1 
0.05 


-0.05 
-0.1 

-0.15 
-0.2 



K,=0.1555, K5e„=0.154 



No Subtraction : 
^ Siihlrartinns i 



500 Configurations and 400 Noise| 



Summation of Insertion Time: (source time+1) to (sink time-1) 
Slope after 3 subtractions (One Source)= 0.0063 ± 0.0089 



10 12 
sinii time (12) 



(c) 



K=0.15S, K,,„=0.154 




No Subtraction ' 
3 Subtractions i 



500 Configurations and 400 Noise] 



Summation of Insertion Time: (source time+1) to (sink time-1) 



Slope after 3 subtractions (One Source)= 0.0034 ± 0.0063 



10 12 
sink time (t2) 



(b) 



I Chirai extrapolation for <x">g (one Nucleon Source) | 



0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 



(d) 



FIG. 11: The ratio of the three-point to two-point functions (summed over insertion time) of the Om — 
\ {Oijj + 04kk) operator, for strange quarks, is plotted against the nucleon sink time t2 at (a) Ky = 0.154 
and Ks = 0.154, (b) Ky — 0.155 and Kg = 0.154, and (c) Ky = 0.1555 and Kg = 0.154. (d) is the hnear 
extrapolation to the chirai limit for the second moment, (a;^) (D.I.), with fixed strange quark plotted against 
(to + ms)a. 



Unfortunately, we did not see any clear signal either for {x'^)u,d (= {x^)u-u = {x^)d-d) 
{x^)s-s, even with unbiased subtractions and multiples sources. The statistical errors are large and 
the error bars overlap either with zero or the signal is at best la to l.ba away from zero for various 
nucleon sources. For the best case (16 nucleon sources), the range for {x'^)u,d is [0.0035,0.0127] and 
that for {x'^)s-s is [0.0016,0.0102]. We conclude that they are consistent with zero in the present 
calculation. 



B. Connected Insertions 



In this section, we will present the results for the first and second moments for the case of 
connected insertions. We will consider both up and down quark currents. Unlike in D.L, they are 
different in the case of C.L As stated earlier, we consider nucleon momentum in the x-direction 
only and we fix the nucleon sink time at ^2 = 16. 
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0.02 



<"» em «■» em ^ 9^- 



-0.02 



-0.04 



-0.06 



-0.08 



k:^=0.154, k^^^=0.154 



1 source 
4 sources ^ ■ 
16 sources ' — e- 




500 Configurations and 400 Noise 




Summation of Insertion Time: (source time+1) to (sink time-1) 


Slope (16 sources)= 0.0055 ± 0.0030 





10 12 
sink time (t2) 



14 



16 



18 



FIG. 12: The ratio of the three-point to two-point functions for the Oiu — ^ {^ijj + ^^kk) operator is plotted 
against the nucleon sink time t2 at k^, = 0.154 and Ks ~ 0.154 for 1, 4, and 16 sources after 3 subtractions. 



0.035 I — 
0.03 
0.025 

0.02 - 
0.015 .. 
0.01 
0.005 - 

— 
-0.005 



Chiral extrapolation for <x 



one source 
4 sources 

16 sources 

fitting line: 1 source 
fitting line: 4 sources 
fitting line: 16 sources 



1" 



(a) 



-0.01 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 
a 



I Chiral extrapolation for <x">g 



one source 
4 sources 
16 sources 
tltting line: 1 source 
fitting line: 4 sources 
fitting line: 16 sources 



0.07 0.08 0.09 



0.11 0.12 0.13 0.14 



(b) 



FIG. 13: Linear extrapolation to the chiral limit for the second moment (a) for {x'^)u,d, and (b) for {x'^)s 
with 1, 4, and 16 nucleon sources. 



First Moments 



First, we will discuss the results for the first moments by using the current O41 for up quarks. 
In Figs. 14(a), 14(b) and 14(c), we plot the ratio in Eq. (22) against the current insertion time ti 
for Ky = 0.154, 0.155, and 0.1555, respectively. In these figures, we see that there is a plateau region 
from the time slice 9 to 13. To extract the values of {x)u+u (CI.) at each k^, we fit a constant 
between the time slice 9 to 13. These values and the corresponding errors are listed in Table X. 

After obtaining {x)u+u (CI.) at finite quark mass, we then linearly extrapolate the valence 
quarks to the chiral limit. Before extrapolation, we have converted the values of {x)u+u (CI.) to 
those of tadpole improved values by using the factors in Eq. (46). As in the case of {x)^^ (D-I-)) 
have extrapolated with the form A + Bruqa. In Table X, we have listed the renormalized linearly 
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1 source 

^INU OUL) j 


1 source 


4 sources 


16 sources 


(Vi T 


Ky — l^sea — 0.154 


0.0020 ± 0.0060 


0.0004 ± 0.0035 


0.0018 ± 0.0021 


0.015 ± 0.0011 


Ky = Ksea ^ 0.155 


0.0056 ± 0.0098 


0.0043 ± 0.0620 


0.0049± 0.0036 


0.0031 ± 0.0020 


i^v — f^sea — U.LOOO 




u.uuyo m u.uuoy 


n 0009 + 00^9 


n 00^9 + 0098 


Linear Extrapolation 


0.0195 ± 0.0213 


0.0168 ± 0.0144 


0.0150 ± 0.0084 


0.0081 ± 0.0046 


{x )s~s 

(D.I.) 


Ky = 0.154, Ksea = 0.154 


0.0020 ± 0.0060 


0.0004 ± 0.0035 


0.0018 ± 0.0021 


0.015 ± 0.0011 


Ky — 0.155, Ksea — 0.154 


0.0048 ± 0.0099 


0.0033 ± 0.0061 


0.0043 ± 0.0036 


0.010 ± 0.0071 


Ky = 0.1555, Ksea = 0.154 


0.0086 ± 0.0134 


0.0061 ± 0.0086 


0.0071 ± 0.0050 


0.0039 ± 0.0027 


Linear Extrapolation 


0.0132 ± 0.0210 


0.0108 ± 0.0137 


0.0113 ± 0.0079 


0.0059 ± 0.0043 



TABLE IX: (x^) (D.I.) for up (down) and strange quarks at various k's and the linearly extrapolated results 
to the chiral limit with different number of nucleon sources. 




10 12 14 

insertion time (t^) 



0.7 
0.6 
0.5 
0.4 
0.3 

0.1 



-0.1 

-0.2 



* * i I z * 



[sink time=16 | 



10 12 14 

insertion time (t^) 



(a) 



(b) 





|K„=0.1555| 




[ 




\ 


sink time=16 | 










{ 



10 12 14 

insertion time (t^) 



0.6 
0.55 

0.5 
0.45 

0.4 
0.35 

0.3 



Ctiiral Extrapolation for <x>y (C.l.) 



(c) 



(d) 



FIG. 14: The ratio of the three-point to two-point functions (C.L), with fixed sink time, <2 = 16, for the ©41 
operator, for up quarks, is plotted against the current insertion time (t) at (a) k„ — 0.154, (b) k„ = 0.155, 
and (c) Ky = 0.1555. (d) is a linear extrapolation to the chiral limit plotted against m^a. 



extrapolated value of {x)u+u (C.L). As stated in Eq. (45), the renormalization factor for this 
operator is 0.972. 

Next, we consider the current for the down quark. Similar procedure is followed as for the up 
quarks to obtain the constants, which give the values of (x)^_(_j's. Again, the fitting is performed 
from the insertion time 9 to 13. The values of {x)j^_^js and their corresponding errors are listed in 
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Table X. 




10 12 14 

insertion time (t^) 




10 12 14 

insertion time (t^) 



(a) 



(b) 



[sink time=16 | 



10 12 14 

insertion time (tj) 



(c) 



0.3 
0.28 
0.26 
0.24 
0.22 

0.2 
0.18 
0.16 
0.14 
0.12 

0.1 



Ciiiral Extrapolation for <x>j (C.I.) 



0.025 0.05 0.075 

m„a 



(d) 



FIG. 15: The ratio of the three-point to two-point functions (C.I.), with fixed sink time, t2 — 16, for the ©41 
operator, for down quarks, is plotted against the current insertion time (t) at (a) Ky — 0.154, (b) Ky — 0.155, 
and (c) Ky — 0.1555. (d) is a linear extrapolation to the chiral limit plotted against m^a. 



Now, we will consider the current O44 — 3 (On -|- O22 + O33) for both up and down quarks. 
The fitting is performed from the insertion time 9 to 11 for both up and down quark currents. 
The values of the {x)u+u and {x)^_^j (CI.) for this operator, along with their errors, are listed in 
Table X. The renormalization factor for this current is 0.953 as obtained from Eq. (45). 



2. Second Moment 



For the second moment, we consider the operator 0411 — ^ (O422 + 0433) for both the up and 
down quarks. We did the similar analysis as in the case of the first moments. The fitting is performed 
from the insertion time 9 to 11 for both up and down quark currents. The values of the {x^)u-u and 
{x'^)d_^ (CI.) for this operator, along with their errors, are listed in Table XI. The renormalization 
factor for this current is 1.116 which we have obtained from Eq. (45). 

We list our results in comparison with those from previous calculations in Table XII. Except 
for our calculations, other groups have averaged the results of the first moments obtained from the 
two operators On and O44 — i (On -|- O22 + O33). As we can see, all the results agree with each 
other within errors. 
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sink time=16 | 
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insertion time (tj) 






(a) 






i i 



0.6 
0.55 
0.5 



~ 0.45 



10 12 14 16 18 20 

insertion time (t.) 



0.4 
0.35 
0.3 



(c) 




10 12 14 16 18 20 

insertion time (t^) 



(b) 



Ciiiral Extrapolation for <x>y (C.I.) 



J !■ 



0.075 



(d) 



FIG. 16: The ratio of the three-point to two-point functions (C.I.), with fixed sink time, t2 — 16, for the 
©44 — i (Oil + O22 + C33) operator, for up quarks, is plotted against the current insertion time ti at (a) 
Ky — 0.154, (b) Ky = 0.155, and (c) k„ = 0.1555. (d) is a hnear extrapolation to the chiral limit plotted 
against m^a. 



For the QCDSF calculations at /x = 5GeV^ (second column), the reported values of (x) are the 
averages of two different procedures in [63, 83], whereas, the result of (x^) is obtained from [83]. 
For both procedures, two different lattices (16^ x 32 and 24^ x 32 ) are used for several k values. For 
the 16^ X 32 lattice, the number of gauge configurations involved is 0(1000) and for the 24'^ x 32 
lattice, the number of gauge configurations involved is O(IOO). All the calculations are performed 
for /3 = 6.0. For the calculations at /i = 4 GeV^ [84] (third column), three k values, 0.155, 0.153, 
and 0.1515, are used on a 16^ x 32 lattice. The number of independent gauge configurations involved 
are 100, 600, and 400 at these n values respectively. All the calculations are performed for (3 = 6.0 
by using standard Wilson action in quenched approximation. 

For the LHPC calculations at = 4 GeV^ (fourth column) [85], a 16'^ x 32 lattice is used for 
three k values. The number of gauge configurations used are 200 at jS = 6.0 for each of the three 
k's in quenched approximation. LHPC calculations are also performed in full (unquenched) QCD 
on a 16^ X 32 lattice. Four different k's are used at (3 = 5.6 for 200 SESAM configurations and 
three different k's are used at (3 = 6.0 for 100 SCRI configurations. 

From this comparison, we see that our results are comparable with other lattice calculations. 
But all the lattice calculations for {x)u-d which involves only the C.I. seem to be larger than the 
experimental result. This is a well known problem and is presumably due to the fact that quark 
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I * * k,=0.154 



-S i 



10 12 14 

insertion time (t.) 



(a) 





10 12 14 16 18 20 

insertion time (t^) 



(b) 



1 1 



* i I * * 



f { 



k,=0.15S5 




10 12 14 

insertion time (t.) 



(c) 



16 18 20 



0.3 
0.28 
0.26 
0.24 
0.22 

0.2 
0.18 
0.16 
0.14 
0.12 

0.1 



Ciiiral Extrapolation for <x>j (C.I.) 



0.025 0.05 
m„ a 



(d) 



0.075 



FIG. 17: The ratio of the three-point to two-point functions (C.I.), with fixed sink time, <2 = 16, for the 
C44 — 5 (Cii + C'22 + ^33) operator, for down quarks, is plotted against the current insertion time ti at (a) 
Ky — 0.154, (b) Ky = 0.155, and (c) k„ = 0.1555. (d) is a Hnear extrapolation to the chiral limit plotted 
against m^a. 



masses are still too heavy compared to the physical ones [79, 80]. 
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{x)t\u (C.I.) 


Ky = 0.154 


0.463 ± 0.014 


Ky = 0.155 


0.448 ± 0.021 


= 0.1555 


0.432 ± 0.027 


Linear Extrapolation 


0.408 ± 0.038 




Ky = 0.154 


0.193 ± 0.007 


Ky = 0.155 


0.179 ± 0.011 


Ky = 0.1555 


0.168 ± 0.014 


Linear Extrapolation 


0.148 ± 0.019 


{x)tU (C.I.) 


Ky = 0.154 


0.491 ± 0.013 


K, = 0.155 


0. IGG ± 0.019 


Ky = 0.1555 


0.458 ± 0.025 


Linear Extrapolation 


0.420 ± 0.035 




Ky = 0.154 


0.224 ± 0.007 


Ky = 0.155 


0.209 ± 0.010 


Ky = 0.1555 


0.203 ± 0.013 


Linear Extrapolation 


0.181 ± 0.018 



TABLE X: Table for the values of {x) (C.I.) for up and down quarks at various kappa values and after linear 
extrapolation to the chiral limit for the O41 and O44 — | (On + O22 + £^33) operators. 



{^')tl\ (CI.) 


Ky = 0.154 


0.128 ± 0.007 


Ky = 0.155 


0.124 ± 0.010 


K^, = 0.1555 


0.122 ± 0.013 


Linear Extrapolation 


0.117 ± 0.018 


i^X-^ (CI-) 


Ky = 0.154 


0.0504 ± 0.0035 


Ky = 0.155 


0.0500 ± 0.0050 


Kt, = 0.1555 


0.0532 ± 0.0066 


Linear Extrapolation 


0.0521 ± 0.0091 



TABLE XI: Tabic for the values of (x^)^^^ (C.I.) for up and down quarks at various kappa values and after 
linear extrapolation to the chiral limit for the O411 — | (C422 + C'433) operator. 
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FIG. 18: The ratio of the three-point to two-point functions (C.I.), with fixed sink time, t2 — 16, for the 
C411 ~ ^ (C422 + C433) operator, for up quarks, is plotted against the current insertion time (t) at (a) 
Ky = 0.154, (b) Ky = 0.155, and (c) k« = 0.1555. (d) is a Hnear extrapolation to the chiral limit plotted 
against mqa. 
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FIG. 19: The ratio of the three-point to two-point functions (C.L), with fixed sink time, t2 = 16, for the 
C411 ~ 5 (C422 + C433) operator, for down quarks, is plotted against the current insertion time ti at (a) 
Ky = 0.154, (b) Ky = 0.155, and (c) k« = 0.1555. (d) is a Hnear extrapolation to the chiral limit plotted 
against m^a. 
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JVEoments 


Kentucky 
(quenched) 
(/x2 = 4GeV2) 


QCDSF 

(quenched) 
{li^ ~ 5GeV2) 


QCDSF 

(quenched) 
{l? = 4GeV2) 


LHPC 

(quenched) 
(/x2 = 4GeV2) 


LHPC 

(/x2 = 4GeV2) 


Experiment 
(m^ = 4GeV2) 




0.408 (38) 
0.420 (35) 


0.410(34) 


0.452 (26) 


0.454 (29) 


0.459 (29) 


0.284 




0.148 (19) 
0.181 (18) 


0.180 (16) 


0.189 (12) 


0.203 (14) 


0.190 (17) 


0.102 


/™2\411 


0.117 (18) 


0.108 (16) 


0.104 (20) 


0.119 (61) 


0.176 (63) 


0.083 




0.052 (9) 


0.036 (8) 


0.037 (10) 


0.029 (32) 


0.031 (30) 


0.025 



TABLE XII: Comparison of results for renormalized first and second moments (connected insertion) in the chiral limit with other lattice calculations and 
phenomenology in MS scheme. 



VII. CONCLUSION AND DISCUSSION 



We have calculated the first and second moments of the proton's parton distribution functions 
for both the connected and disconnected insertions in lattice QCD. The lattice calculations are 
carried out on quenched 16^ x 24 lattices with (3 = 6.0 and quark masses which correspond to pion 
masses of 650(3), 538(4) and 478(4) MeV, and nucleon masses at 1291(9), 1159(11) and 1093(13) 
MeV respectively. The physical results are obtained from linear extrapolation to the physical point 
(to the chiral limit). The connected insertion results turn out to be consistent with the previous 
quenched and full QCD calculations. 

The calculation of the moments for the sea quark distribution in the D.I. is carried out for the 
first time. With 400 Z2 noise, 500 gauge configurations, and with unbiased subtractions and 16 
nucleon sources, we are able to obtain results with ~ 5(T signals for the first moments. The result 
of {x)s = 0.027 lb 0.006 can be used to constrain {x)s in the phenomenological fitting of parton 
distribution functions which is uncertain in the range 0.018 < {x)s < 0.04 [8]. More interestingly, 
we find that 



{x) 



s+s 



{x) u+u 



dx X {s{x) + s(x)) 



D.I. 



+.(.)) 



= 0.88 ± 0.07, (48) 

D.I. 



which is about twice as large as the average phenomenological value from fitting the parton distri- 
bution functions to experiments [8] 

dx X [s{x) + s{x)] 

0.27 - 0.67. (49) 



dx X [n(x) + d{x)] 

This difference is understandable and has been anticipated from the path-integral formulation 
of parton degrees of freedom [21]. The ratio in the lattice calculation involves u/d quarks in the 
disconnected insertion (quark loops) [89], while the phenomenological ratio involves the u{x) + d{x) 
in the connected insertion as well. The discrepancy suggests that the momentum carried by the 
u/d quarks are roughly equally shared in the connected sea and the disconnected sea. The fact 
that there is anti-quark in the connected sea is demonstrated by the large u{x) — d{x) difference 
from the Gottfried sum rule violation [39] which cannot be accommodated by the small u and 
d difference due to isospin breaking in the disconnected sea. The combined Gottfried sum rule 
violation and the discrepancy of the ratios in Eqs. (48) and (49) suggests the following form for 
the anti-parton distribution functions [21]: 

u{x) = aox'^^il-xf^ +hux^^{l-xf^, 
d{x) = aox^'il-x)"^ +bdx^^{l-x)^\ 

s{x) = asx^^il-xf^ (50) 

where the first terms are for the disconnected seas with pomeron exchanges so that ai ~ — 1 and 
the second terms are from the connected sea with reggeon exchanges so that bi ~ — a_R. The 

present lattice calculation suggests that — can be constrained to the ratio in Eq. (48). We should 

ao 

emphasize that the current result is based on a quenched lattice calculation with linear chiral 
extrapolation from relatively heavy quarks. We should take it with a sizable grain of salt. As far 
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as phenomenological fittings are concerned, the u{x) — d{x) has been taken into account which has 
a smah x behavior of ~ 3;~^/^[86]. However, the conventional ansatz s{x) ~ u{x) + d{x) used in 
recent fittings [8, 9] is obviously inadequate [21]. One needs to differentiate the different small x 
behaviors in the connected sea and disconnected sea and fit the and s{x) accordingly. 

We have made an attempt to calculate (x^)^ to see if it is s{x) or s{x) which is leading in 
large x. Our result has a tendency to be positive, similar to the experimental tendency of s_ = 

j dxx {s{x) — s{x)) being positive, but it is consistent with zero within error. We will see if the 

signal is stronger with the inclusion of 2+1 flavor dynamical fermions. 

So far, our results are obtained in the quenched approximation with relatively large quark masses 
and small volume. They are subjected to large systematic errors. We will focus our attention next 
on the dynamical fermion calculation with 2+1 flavor dynamical clover fermion configurations [87], 
and systematically move to smaller quark masses, larger volumes, and the continuum limit [88]. 
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APPENDIX A: THREE-POINT CORRELATION FUNCTIONS 
1. General Considerations 

For the forward matrix elements, the three point-function for any general operator O (color and 
spin indices suppressed) is defined as 

G°^o^(t2,ii,p1 = Y. e-^^^-^"^^~"^"Mo|T(x"(x2)0(xi)x^(xo))|0), (Al) 

where t = f2 is the nucleon sink time, t = ii is the current insertion time, t = to is the nucleon 
source time, and pi is the momentum of the nucleon, respectively. We now consider the following 
situations: 

• When t2 > ti > Iq, we get 

= 5^ e-^^^-(^^^-^^o) E (0 I 1^2,92, S2)(n2,g-2,52|0(xi)|ni,gi,si) 

{nuqi,si\ x^(xo) |0) 
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= N'Yl 6-^°^^*^-*°^ ( 1 x"(xo) I 0+,p, s){0+,p,s\ O{xo) I 0+,p, s) 

s 

{0+,p,s\x^{xo)\0) 
+ iV2^e-^p"(*^-*°)(0|x"(a:o)|0-,p,s)(0-,p,s|0(xo)|0-,p,s) 

s 

{0-,p,s\x^{xo)\0)+ J2 [e-<'"(*^-*i)e-<'(*i-*o)c(^)W"^(p-) 

e=+,- ^ 



(A2) 



where N is the number of lattice sites, and the superscript + (— ) represents positive (negative) 
parity state. 

CW(i)a/3(^-^ ^ N^{0\x''{xo)\0^'\p,s){o('\p,s\0{xo)\l^'\p,s) 
(lW,p,s|x^(a;o)|0), 

Cim)aP^^ ^ N^{0\^-{xo)\l'^'\p,s){l^'\p,s\0{xo)\o'^'\p,s) 
(0W,p,s|x^(xo)|0), 

/(i)°/^(nf ) , nP , p) = iV' 5^ ( I (xo) I nf ^ ,p,s){ , p, s | O(xo) | nf ^ , p, s) 

s 

{nf\p,s\x''{xo)\0). (A3) 

• When ti > t2 > to, we get 

ni,n2,s 

(n2,0,s|x'*(a;o) {ni,p,s\x^{xo)\0) 
= 0. since, for a twist-two operator, ( | O{xo) | n2, 0, s ) =0 (A4) 

• When t2 > to > ti, we get 

G%%^{t2,ti,p) = Yl e-^^"'^*^-*°^e-<'^*°-*i)(0|x"(xo)|n2,p,s) 

ni,n2,s 



(n2,P,s|x'^(a;o) |ni,0,s) (ni,0,s|0(a;o) |0) 



= 0. 



since, for a twist-two operator, ( ni, 0, s | O{xo) | ) = (A5) 
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2. Ratios for Disconnected Insertion 



We will show a sample calculation for extracting the disconnected first and second moments by 
considering the O = On operator, where i = 1,2,3. The calculation for the other operators will 
be similar. For this operator, we will use the nucleon with one unit of lattice momentum in the 
i-ih direction. First consider the first and the second terms of Eq. (A2) only. Taking the trace with 
parity projection operator, F^" in Eq. (A2), we get 

F^" X [first term + second term of Eq. (A2)] 



^2 ^g-EO+(t2-to) 



2 a I m 



(2k) i VA^^r; ' (2k) i \NE^p 



^+|2 (t2-*o) (i + a^ \ Mzi 

(2k)3 ' E^- m+ I 2k 



(A6) 



where (x)"*" = {x) is the first moment, (x) is an unknown constant. 

We will now sum Tr \^ GNON{'t2,ti^'P)\ over the current insertion time, ti, from an initial time, 
= ^0 + 1) to a final time, tj = ^2 ~ 1) so that, tf > ti > to, where the nucleon source is at to, and 
t2 is the sink time. Then, by using Eqs. (A2), (A3) and (A6), we get 

*/ */ 

^ Tr [TGNoMt2,h,pO] = ^F^-G^^^^^^(t2,ti,p1 



a m 



■Ti. (2^)' E', 



|^|2g-i?0(t2-to) 



1 + 



m 



X) Pi 



Et~ m / 2k 



+ E 



+ e" 



^^'*'(t2-ti)g-£;r'(ti-t«)c.{W2)(^ 



-,0(6) 



OO r (6') (9) 



(A7) 
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Dividing by the two-point function, we get 



where 



E 

t\=ti 



Tr [GNOuN{t2,ti,p)] 



{x)pi ^ _ {x)pi 
2k ^ 2k 



d'=+-n'(0')=l 



+ 



e=+ -n(9)=l 



m=l 



ki 



(2^)3 E^^ 

k2 = E 

k3= E 



constant), 



1 



c, 



{em, 

ii 



_ 1 — Q-\l^v ~'^P ) 



constant), 



E 



C7if)(2)(= constant). 



1(9) 0(6)^ ^ii'"' "*(= constant). 



l_e(<'-<^ ) 



/«(nW,nP,p-) 



From Eq. (A8), we get 



V- Tr [GNoMt2,ti,p)] - to) » 1 {x)pi 

y, rr^ r-r^ ^ T, ~ ^ ~ l2+COnSt., 



tl=ti 



Tr[TGNNit2 



2k 
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where 



const. = — 



{x)Pi 
2k 



ti + kik2 - kik5. 



(All) 



APPENDIX B: UNBIASED SUBTRACTION 

The set of traceless matrices is obtained from the hopping parameter expansion of the inverse 
of the fermion matrix M which is given by 



where 



M-^ = I + kD + k^D^ + K^D^ + • • • , 



Dx,y = I] (1 - 7M)a/3 Uf{x) d^,y-a, + (1 + 7n)aP C/J^^X - a^) d^,y+a^ 



(Bl) 



(B2) 



X, y arc space-time co-ordinates, a, (3 are spin indices, and a, b are color indices. If for some currents, 
some powers of D are not traceless, we make them traceless by subtracting their traces from 
themselves. Now, we are going to calculate the traces for all the three currents. We suppress the 
argument U from each of the propagators for convenience. 



1. Two- Index Operators 



The disconnected part of O^i, is given by 
Loop 

= + - 1V[5W'""(xi-a,,xi)7^C/^'"(xi-a,)] 



XI 



+ Tr [-S(^)-"(a;i, XI - a,) 7^ Uri^i - - 1^ [-S(^)"^"(a;i, xi + a,) 7^ ^^^^'"(^i)] 



(B3) 



where A = — i for = 4, v = 1,2,3; X = +1 for = = 1, 2, 3 and A = — 1 for = 1/ = 4. Let us 
consider the first term of the loop 



8a 



^ TV \s(f^^^{xi + a,, xi) 7^ Urixi) 



'^12 Yl ('^li,n,r Xxi+a.,m,^) {7ti)r^ K^'ixi) [using noise] 



XI m,£, 



^ E E 7 E ^'ln,r mI-'^-\x, + a., .) ry^,,,, (7,)., un^,). 



x\ m,^,k,cr 1=1 



(B4) 
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In Eq. (B4), we will substitute the hopping parameter expansion of M ^. Then, from Eq. (B4), 
the term becomes 



1 ^ 



L 

x\ 1=1 



(B5) 



Similarly, if we consider the fourth term of the Loop in Eq. (B3), we get 

A 



8a 



8a L 

Xl (=1 



(B6) 



Combining Eqs. (B5) and (B6), we find that 

A r, 



8a 



Tr + a„xi)7M^r(^i)] " Tr [^(^^'"^xi, xi + a,) 7^ C/t"-(xi)] 



^E E 7 E 2Im kt,^^,,,^,^^^^„^^^(^^)^^C/-(xO' 



8a ^ ^ L 

XI m,§ (=1 



(B7) 



By considering similar combinations of other terms, we can prove that for the O^^ operator, the 
real part of the loop is zero for the first term, I, and for the Oi^i operator, the imaginary part is 
zero. So, we cannot consider the first term, I. If we substitute M~'^ by k?D^ and k'^D^, the trace 
is going to be zero, since there will be no plaquette term. But if we substitute by kD and 

K^D^, there is a possibility of having a non-zero trace. Since the operators, and O^^, can have 
different traces (or, no traces), we will consider them separately. 

• First term of 044: 



E [5(^)"^'^(xi + ai, Xl) 74 Urixi) 



8a 

—i 
8a 



Xl 



E E 7E^iI,n,.M€t'^"''(^i + «.,^)<fc,.(74).cc/r(^i)- 



L 

XI m,^,k,a 1=1 



(B8) 



If we substitute M ^(xi + at, z) in Eq. (B8) by kD, the only non- vanishing term could be. 



X ^fi) ^x,y+afj 



So we get 



^E^r [5(^)"^"(xi + ai,xi)74C/r(xi)' 



XI 
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-i 1 

xi m,^,k,(j 1=1 



11=1 



8a 



0. 



'^(n,fc)(T,f7) (74)r5(l + 7i)C(TC/'/"''(a:i)f/"'"(xi) ( putting // = z ) 



(B9) 



We now replace by k^D^, the non-vanishing terms could be 

^ {(l-7M')(l + 7^')(l + 7p)}a/3{C^/.'(a;)C/^,(x + a^/ -a^O 

+ {(l + 7p')(l-7i^')(l+7p)}a/3{f^^/(2;- v)^'^'(^~ V)^p(^~ V +«'^' -«p)} 



{(1 + 7/,')(l +7i'')(l - 7p)}a/3{f^^'(a; - ai,')ul,{x - a^' - a,,/)f7p(x - a^/ - a,,/)}" 



Jx-a i-a^i ,y-ap 



(BIO) 



Let us consider the first term of Eq. (BIO). 

-^^Tr [5(^)'""(xi+a„xi)74C/r(xi)' 



8a 

—i 
8a' 



XI 



1 ^ 



L 

a;i m,i;,k,cr 1=1 



{(1 - 7^0(1 + 7i.')(l + 7p)ka{f/M'(2;i + «i) 



Ul,{xi + ai + a^' - Oi,') Uj,{xi + + a^' - Oj,/ - Cp)}™'' J^-i+ai+a ,-a^,,z+ap 



<fc,.(74).c^r(xi)- 



(Bll) 



In order to have non-vanishing delta function, we must have either /i' = p, v' = i; or fi' = v', p = i; 
or p' = p = u' = i. Let's take the case p' = p,v' = i. Then the first term of Eq. (BIO) gives 



''^^Yl Y IZ^(".fc)(^.'^)(^4)rc(l-7p')(l+7i)(l+7p')5<T{f^p'(a;i + 



ai 



xi m,^,k,cr fi' 

U}{x, + a,,)Ul,{x,)r'ur{xi) 
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8a ^ 
0. 



(B12) 



For other two possibilities fi' = v', p = i and n' = p = v' = i, the trace for the spin part is zero 
due to the multiphcation of 7 matrices. 

It can be shown that second and third terms of Eq. (BIO) also have zero traces. In a similar 
manner, we can show that other terms of the operator O^i also have zero trace. 

Now, we will consider the operator, O^/^. 

• First term of Of^^: 



-J^Tr [5(/)™"(xi + a^,xi)7Mf/r(^i) 



A 



_E E iE^iLn,.M(;^)™^xi + a^,z)<,,,(7^).5[/;-(xi), (B13) 

xi m,t^,k,a 1=1 

where A = +1 for p = v = 1,2,2> and A = —1 for p = u = A. 

If we substitute M^^{xi + a^, z) in Eq. (B13) by kD, the only non-vanishing term could be 

4 



So we get 



8a 



^Tr [5(^)™"(xi+a^,xi)7Mf/; 



Xl 



8a ""E E ^ E ^-L-,- 

Xl m,^,k,a 1=1 

vUA7,Uu;^{xi) 



p=i 



= ^'^Yl E ^(n,fc)(ra)(7/.)r5(l + 7Mk-^M"'^(a^l)f^M"'(^l) ( puttiug = ^ ) 

Xl m,^,k,a 

8a 

And, if we substitute M~^(xi + a^, z) by k^D^, the non- vanishing terms could be 
E {(l-V)(l+7i^')(l+7p)W{f^M'(3;)^^^/(x + a^/-a^/) 



(B14) 
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{(1 + 7m')(1 - 7i^')(l + lp)}ai3{Ul,{x - ap>)Uyi{x - af,')Uj,{x - a^/ + a,,' - Op)}" 

^x—ai+a^,,y+ap 



{(1 + 7/»')(l +7i^')(l - 7p)}af3{Ul,{x - ap')ul,{x - a^' - ayi)Up{x - a^/ - a^')}" 



Jx-a ,y-ap 



(B15) 



Let us consider the first term of Eq. (B15) 
A 



8a 



8a ^ ^ L 

x\ m,^,k,a 1=1 



{(1 - 7pO(l + 7v')(l + lp)}^a{Ut.'{xi + a^) 



Ul,{xi + a^ + a^/ - a^/) ?7T(xi + + a^/ - a,,/ - Op)}™ (5xi+a^+a^,-a,y ,2+ap 
<.,.(7.).cC/r(^i)- 



(B16) 



In order to have non- vanishing delta function, we must have either fj,' = p, u' = fi, or ji' = v' , p = ^, 
or p' = p = u' = p. Let's take the case p' = p,v' = p. Then the first term of Eq. (B15) gives 



'^^Yl Y Y^{n,k){T,a){lt,)Ti{'^-lp'){'^ + lt,){'^+lp')ia{U^,'{,Xi+a^) 



X_ 

c/t(xi+a^o^^'(^i)r'c^r(^i) 

XI m n'^fj. 



(B17) 



For other two possibihties p' = v', p = p and p' = p = u' = p, the trace for the spin part is zero 
due to the muhiphcation of 7 matrices. Let now us consider the second term of Eq. (B15). Then 



A 

8a 



^ Tr + ap, xi) 7^ UJi'^ixi) 



^ ^«^'E E rE^itn,^ E + 7.0(1- 7.0(1 +7p)Mf^^,(^i+ a, -a,0 



8a'" ^ ^ L 

Xi m,i;,k,cr 1=1 



Uy'{xi + ap - a^/) C/^(xi + a^ - a^/ + a^,/ - ap)}"'^4i+a^-a^,+a„,,z+ap 



<.,.(7Mkf/r(^i)- 



(B18) 
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In order to have non-vanishing delta function, we must have either /j,' = fi, p = v', or /i' = i'', p = p,, 
or p' = p = u' = p. Let's take the case p' = p, v' = p. Then the second term of Eq. (B15) gives 



A 

8a 



(7M)r5 {(1 + 7m)(1 - 7.0(1 + lu')]ia{Ul{xi) U,,{xi) 

Xi m,i;,k,cr fl' 



0. 



(B19) 



For other two possibihties p' = v', p = p, and p' = p = v' = p^ the traces for the spin part are 
zero due to the multiphcation of 7 matrices. Let us consider the third term of Eq. (B15). Then 



8a 

X_ 
8a' 



IV [^(^^'""(xi + a^, xi) 7^ U;^{xi) 



1 ^ 



afi m,^,k,<j 1=1 



{(1 + 7^0(1 + 7v')(l - 7p)}5a{C/^/(xi + a^ - a^O 



C/J,(xi + a^ - ttfj,' - ayi) Up{xi + a^ - a^i - ay')}'^ 5xi+a^-a^,-a^,,z-a 



(B20) 



In order to have non-vanishing delta function, we must have either p' = p, v' = p, or p' = p, p = 
v', or p' = p = u' = p. Let's take the case p' = p, v' = p. Then the third term of Eq. (B15) gives 



X] X] X]'^(".'^)(^.'^)(>)^«(^+V)(l + 7A.)(l-7M')5'x{f^^'(2;i+a;,-a^/) 



x\ m,t^,k,cr fi' 

X 



8a 



For other two possibilities p' = p, p = v' , and p' = p = v' = p, the traces for the spin part are 
zero due to the multiplication of 7 matrices. 



Doing similar calculations for other terms of O^^, we get the trace for the kD term 

Trace = —2k [l2V^ - {-12V^) + 12V'^ - (-121/3)1 = AgeK^^^ 
8a 8a 

and for the k^D^ term 



(B22) 



Trace ~ ^I^k^E E E 



{C/^/(xi + a^)?7^(xi + a^>)Ul,{xi)Uf,{xi) 



+ Ul,{xi + a^ - a^i) UUxi - a^/)C/^/(a;i - a^/)f7^(xi)} 
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+ {Ui,'{xi - a^)U^{xi - + a/,/)C/^,(xi)[/^(xi - a^) + - a^- a^O 



(B23) 



2. Three-Index Operators 

As in the case of two-index operator, we can show that the first term I of the hopping expansion 
for tliis operator is real. But the loop part is imaginary. So, we can not consider the first term. If 
we substitute by kD and k^D^, the trace is going to be zero, since there will be no plaquette 
terms. But if we substitute by k^D^, there is a possibility of non-zero trace. One of the terms 
in the loop part of Om is 



-1 
24^ 



Tr 



+ a4 + ai,xi) (7^) \ (xi) uf'^ixi + 04) + U^' {xi) 



UT'^ixi + ai 



-1 1 

XI m,^,k,a 1=1 

( U2"''{xi)Ur'"'ixi + 04) + Ur' {xipf^'ixi + ai). (B24) 
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If we substitute mL ^^"^'^ by k^D^ term, the only non- vanishing term could be 



4 4 



fx' = l v' = l ^ 



Therefore, we get 
-1 



24.0? 



+ Unxi)Ur'^{xi+ai)) 



-1 1 ^ ^ 



{(1 + 7^0(1 + 7^')}«a 

{?7^,(xi + a4 + ttj - a^')^r!'(a^l + 04 + - V - a,,')}"*''(5j;i+a4+a;_a^,,2+a^, 



-1 



24a2- E E Wr,a)(7ik{(l + 7i)(l + 74)W{i7/(a;i+a4))C/i(xi)r^ 
{U^ixipiixi + 04)}"'" + {?7i(xi)C/4(xi + ai)}""^! (putting /x' = i and z/' = 4) 



+ E ^n,k)ir,a){l^U{i^+l^)i^+l^)hAUKxl+a,))Uj{x,)^'^ 



{U^ixipiixi + a4)r'" + {Uiixi)U4ixi + ai)y 



(putting /x' = 4 and v' = i) 



8^«'^ 5^Re{{C//(xi + a4)C/i(xi)f/i(xi)C/4(xi + ai)r'" + 3 



XI 



(B25) 



Calculating all the other terms in a similar manner, we can show that 

Trace = ^Sk^ ^ ^ Re {{f//(xi + a4)?7](xi);7i(xi);74(xi + aO}"*"* + 3 



- 8k^Y1 ^lra{ul{xi-ai)Ui{xi-ai)U4{xi)Uj{xi + a4-ai)r 

XI m 

- 8k^ ^Im{f/4(xi - 04 + ai)U^i {xi)ul{xi - a4)Ui{xi - 04)^ 



x\ 
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— Re {Ui{xi — 04 — ai)U4{xi — 04)^/ (xi — aj) 

Ulixi - 04 - ai)r"' + 3 
+ Re {C/i(a;i — 04 — aj)C/4(xi — a4)?7/(xi 

-04-00}"™ + 3 

— ^^lm{U4{xi — + ai)uj {xi)ul{xi — a4)Ui{xi — a^)} 

- 8k^^ ^Im{C/](xi - ai)?7j(xi - ai)?74(xi)C//(xi + 04 - Oj)}'^ 

- 8^2^ ^Re{{uUxi+a^)ul{xi)U,{xi)U^{xi+ai)r"' + 3 
^ ^^^^^^ ^Im{[/](xi - ai)C/i(xi - ai)C/4(xi)C//(xi + a4 - Oi)} 



xi 



+ ^Im{?74(xi - 04 + aj)C//(xi)C/|(xi - a4)?7j(xi - 04)}' 



(B26) 
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